Explore CONUS404 Dataset

This dataset was created by extracting specified variables from a collection of wrf2d output files, rechunking to better facilitate data extraction for a variety of use cases, and adding CF conventions to allow easier analysis, visualization and data extraction using Xarray and Holoviz.

import os
os.environ['USE_PYGEOS'] = '0'

import fsspec
import xarray as xr
import hvplot.xarray
import intake
import metpy
import cartopy.crs as ccrs

1) Open dataset from Intake Catalog

• Select on-prem dataset from /caldera if running on prem (Denali/Tallgrass)

• Select cloud/osn data on S3 if running elsewhere.

**NOTE This notebook reads data from OSN by default, which is free to access from any environment. If you change this notebook to use one of the CONUS404 datasets stored on S3 (options ending in -cloud), you will be pulling data from a requester pays bucket. This means you have to set up your AWS credentials, else we won’t be able to load the data. Please note that reading the -cloud data from S3 may incur charges if you are reading data outside of the us-west-2 region or running the notebook outside of the cloud altogether. If you would like to access one of the -cloud options, add a code cell and run the following code snippet to set up your AWS credentials:

os.environ['AWS_PROFILE'] = 'default'
%run ../environment_set_up/Help_AWS_Credentials.ipynb

# open the hytest data intake catalog
hytest_cat = intake.open_catalog("https://raw.githubusercontent.com/hytest-org/hytest/main/dataset_catalog/hytest_intake_catalog.yml")
list(hytest_cat)

['conus404-drb-eval-tutorial-catalog',
 'nhm-v1.0-daymet-catalog',
 'nhm-v1.1-c404-bc-catalog',
 'nhm-v1.1-gridmet-catalog',
 'conus404-catalog',
 'nwis-streamflow-usgs-gages-onprem',
 'nwis-streamflow-usgs-gages-cloud',
 'nwm21-streamflow-usgs-gages-onprem',
 'nwm21-streamflow-usgs-gages-cloud',
 'nwm21-streamflow-cloud',
 'nwm21-scores',
 'lcmap-cloud',
 'rechunking-tutorial-cloud']

# open the conus404 sub-catalog
cat = hytest_cat['conus404-catalog']
list(cat)

['conus404-hourly-onprem',
 'conus404-hourly-cloud',
 'conus404-hourly-osn',
 'conus404-daily-diagnostic-onprem',
 'conus404-daily-diagnostic-cloud',
 'conus404-daily-onprem',
 'conus404-daily-cloud',
 'conus404-monthly-onprem',
 'conus404-monthly-cloud']

## Select the dataset you want to read into your notebook and preview its metadata
dataset = 'conus404-hourly-osn' 
cat[dataset]

conus404-hourly-osn:
  args:
    consolidated: true
    storage_options:
      anon: true
      client_kwargs:
        endpoint_url: https://renc.osn.xsede.org
      requester_pays: false
    urlpath: s3://rsignellbucket2/hytest/conus404/conus404_hourly_202302.zarr
  description: 'CONUS404 Hydro Variable subset, 40 years of hourly values. These files
    were created wrfout model output files (see ScienceBase data release for more
    details: https://www.sciencebase.gov/catalog/item/6372cd09d34ed907bf6c6ab1). This
    dataset is stored on AWS S3 cloud storage in a requester-pays bucket. You can
    work with this data for free in any environment (there are no egress fees).'
  driver: intake_xarray.xzarr.ZarrSource
  metadata:
    catalog_dir: https://raw.githubusercontent.com/hytest-org/hytest/main/dataset_catalog/subcatalogs

2) Parallelize with Dask

Some of the steps we will take are aware of parallel clustered compute environments using dask. We’re going to start a cluster now so that future steps can take advantage of this ability.

This is an optional step, but speed ups data loading significantly, especially when accessing data from the cloud.

%run ../environment_set_up/Start_Dask_Cluster_Nebari.ipynb
## If this notebook is not being run on Nebari/ESIP, replace the above 
## path name with a helper appropriate to your compute environment.  Examples:
# %run ../environment_set_up/Start_Dask_Cluster_Denali.ipynb
# %run ../environment_set_up/Start_Dask_Cluster_Tallgrass.ipynb

The 'cluster' object can be used to adjust cluster behavior.  i.e. 'cluster.adapt(minimum=10)'
The 'client' object can be used to directly interact with the cluster.  i.e. 'client.submit(func)' 
The link to view the client dashboard is:
>  https://nebari.esipfed.org/gateway/clusters/dev.0ad47654db8042b393bf20edd4b1db84/status

3) Explore the dataset

print(f"Reading {dataset} metadata...", end='')
ds = cat[dataset].to_dask().metpy.parse_cf()
print("done")
# Examine the grid data structure for SNOW: 
ds.SNOW

Reading conus404-hourly-osn metadata...

done

<xarray.DataArray 'SNOW' (time: 368064, y: 1015, x: 1367)>
dask.array<open_dataset-986a06d6b28d681a9e9ca356f92aad3eSNOW, shape=(368064, 1015, 1367), dtype=float32, chunksize=(144, 175, 175), chunktype=numpy.ndarray>
Coordinates:
    lat        (y, x) float32 dask.array<chunksize=(175, 175), meta=np.ndarray>
    lon        (y, x) float32 dask.array<chunksize=(175, 175), meta=np.ndarray>
  * time       (time) datetime64[ns] 1979-10-01 ... 2021-09-25T23:00:00
  * x          (x) float64 -2.732e+06 -2.728e+06 ... 2.728e+06 2.732e+06
  * y          (y) float64 -2.028e+06 -2.024e+06 ... 2.024e+06 2.028e+06
    metpy_crs  object Projection: lambert_conformal_conic
Attributes:
    description:   SNOW WATER EQUIVALENT
    grid_mapping:  crs
    long_name:     Snow water equivalent
    units:         kg m-2

xarray.DataArray

'SNOW'

• time: 368064
• y: 1015
• x: 1367

• dask.array<chunksize=(144, 175, 175), meta=np.ndarray>

  Array Chunk Bytes 1.86 TiB 16.82 MiB Shape (368064, 1015, 1367) (144, 175, 175) Dask graph 122688 chunks in 2 graph layers Data type float32 numpy.ndarray

• Coordinates: (6)
  • lat
    (y, x)
    float32
    dask.array<chunksize=(175, 175), meta=np.ndarray>

    description :

    LATITUDE, SOUTH IS NEGATIVE

    long_name :

    Latitude, south is negative

    standard_name :

    latitude

    units :

    degree_north

    Array Chunk Bytes 5.29 MiB 119.63 kiB Shape (1015, 1367) (175, 175) Dask graph 48 chunks in 329 graph layers Data type float32 numpy.ndarray

  • lon
    (y, x)
    float32
    dask.array<chunksize=(175, 175), meta=np.ndarray>

    description :

    LONGITUDE, WEST IS NEGATIVE

    long_name :

    Longitude, west is negative

    standard_name :

    longitude

    units :

    degree_east

    Array Chunk Bytes 5.29 MiB 119.63 kiB Shape (1015, 1367) (175, 175) Dask graph 48 chunks in 329 graph layers Data type float32 numpy.ndarray

  • time
    (time)
    datetime64[ns]
    1979-10-01 ... 2021-09-25T23:00:00

    standard_name :

    time

    array(['1979-10-01T00:00:00.000000000', '1979-10-01T01:00:00.000000000',
           '1979-10-01T02:00:00.000000000', ..., '2021-09-25T21:00:00.000000000',
           '2021-09-25T22:00:00.000000000', '2021-09-25T23:00:00.000000000'],
          dtype='datetime64[ns]')

  • x
    (x)
    float64
    -2.732e+06 -2.728e+06 ... 2.732e+06

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    meters

    array([-2732097.901154, -2728097.901154, -2724097.901154, ...,  2723902.098846,
            2727902.098846,  2731902.098846])

  • y
    (y)
    float64
    -2.028e+06 -2.024e+06 ... 2.028e+06

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    meters

    array([-2027960.899637, -2023960.899637, -2019960.899637, ...,  2020039.100363,
            2024039.100363,  2028039.100363])

  • metpy_crs
    ()
    object
    Projection: lambert_conformal_conic

    array(<metpy.plots.mapping.CFProjection object at 0x7f3295dc1f90>,
          dtype=object)

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1979-10-01 00:00:00', '1979-10-01 01:00:00',
                   '1979-10-01 02:00:00', '1979-10-01 03:00:00',
                   '1979-10-01 04:00:00', '1979-10-01 05:00:00',
                   '1979-10-01 06:00:00', '1979-10-01 07:00:00',
                   '1979-10-01 08:00:00', '1979-10-01 09:00:00',
                   ...
                   '2021-09-25 14:00:00', '2021-09-25 15:00:00',
                   '2021-09-25 16:00:00', '2021-09-25 17:00:00',
                   '2021-09-25 18:00:00', '2021-09-25 19:00:00',
                   '2021-09-25 20:00:00', '2021-09-25 21:00:00',
                   '2021-09-25 22:00:00', '2021-09-25 23:00:00'],
                  dtype='datetime64[ns]', name='time', length=368064, freq=None))

  • x
    PandasIndex

    PandasIndex(Float64Index([-2732097.901153542, -2728097.901153542, -2724097.901153542,
                  -2720097.901153542, -2716097.901153542, -2712097.901153542,
                  -2708097.901153542, -2704097.901153542, -2700097.901153542,
                  -2696097.901153542,
                  ...
                   2695902.098846458,  2699902.098846458,  2703902.098846458,
                   2707902.098846458,  2711902.098846458,  2715902.098846458,
                   2719902.098846458,  2723902.098846458,  2727902.098846458,
                   2731902.098846458],
                 dtype='float64', name='x', length=1367))

  • y
    PandasIndex

    PandasIndex(Float64Index([-2027960.8996368449, -2023960.8996368449, -2019960.8996368449,
                  -2015960.8996368449, -2011960.8996368449, -2007960.8996368449,
                  -2003960.8996368449, -1999960.8996368449, -1995960.8996368449,
                  -1991960.8996368449,
                  ...
                   1992039.1003631551,  1996039.1003631551,  2000039.1003631551,
                   2004039.1003631551,  2008039.1003631551,  2012039.1003631551,
                   2016039.1003631551,  2020039.1003631551,  2024039.1003631551,
                   2028039.1003631551],
                 dtype='float64', name='y', length=1015))

• Attributes: (4)

  description :

  SNOW WATER EQUIVALENT

  grid_mapping :

  crs

  long_name :

  Snow water equivalent

  units :

  kg m-2

Looks like this dataset is organized in three coordinates (x, y, and time). There is a metpy_crs attached:

crs = ds['SNOW'].metpy.cartopy_crs
crs

<cartopy.crs.LambertConformal object at 0x7f3295a0cd00>

Example A: Load the entire spatial domain for a variable at a specific time step

%%time
da = ds.SNOW_ACC_NC.sel(time='2009-12-24 00:00').load()
### NOTE: the `load()` is dask-aware, so will operate in parallel if
### a cluster has been started. 

CPU times: user 770 ms, sys: 216 ms, total: 986 ms
Wall time: 20 s

da.hvplot.quadmesh(x='lon', y='lat', rasterize=True, geo=True, tiles='OSM', cmap='viridis').opts('Image', alpha=0.5)

Example B: Load a time series for a variable at a specific grid cell for a specified time range

SIDE NOTE To identify a point, we will start with its lat/lon coordinates. But the data is in Lambert Conformal Conic… need to re-project/transform using the built-in crs we examined earlier:

lat,lon = 39.978322,-105.2772194    
x, y = crs.transform_point(lon, lat, src_crs=ccrs.PlateCarree())   
print(x,y) # these vals are in LCC

-618215.7570892666 121899.89692719541

%%time
da = ds.PREC_ACC_NC.sel(x=x, y=y, method='nearest').sel(time=slice('2013-01-01 00:00','2013-12-31 00:00')).load()

CPU times: user 208 ms, sys: 20.5 ms, total: 229 ms
Wall time: 20.8 s

da.hvplot(x='time', grid=True)

Stop cluster

client.close(); cluster.shutdown()

/home/conda/users/ef214d1986de9477d5223419ba4838d76e9839b19e4fb02c9e208c97cab5a61b-20230329-131713-443303-165-pangeofu/lib/python3.10/site-packages/dask_gateway/client.py:1014: RuntimeWarning: coroutine 'rpc.close_rpc' was never awaited
  self.scheduler_comm.close_rpc()